Programmer Power Tools/Programmer Power Tools.iso/microcrn/issue_40.arc/DAIMS.ARC/MATRIX1.CPP

C/C++ Source or Header | 1988-02-10 | 14.6 KB | 630 lines

#include "matrix.hxx"
#include <math.h>
#include <ctype.h>
#include <time.h>
#ifndef DOS
#include <sys/param.h>
#endif
#define UPCASE(X)  (X = islower(X) ? toupper(X) : X) /* force a char to upper case */

/*
-*++ matrix::matrix(): contructs a matrix of any size and inits to 0
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

  matrix::matrix(int rows = 1 , int cols = 1, double initval = 0)
{
  p = new matrep;
  p->m_vec = new _Vector*[rows];
  p->m_rows = rows;
  p->m_cols = cols;
  for (int i = 0; i < rows; i++)
    p->m_vec[i] = new _Vector(cols);
  for (i=0; i < rows; i++)
    for (int j = 0; j < cols; j++)
      (*this)[i][j] = initval;
  p->n = 1;            /* so far, only one reference to this data */
}


/*
-*++ matrix::matrix(): adds a new reference to an existing matrix
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed: "Reference counting" is a method of managing garbage.
** Instead of creating a whole new object every time you use (for instance) an
** equal sign, it makes more sense just to point at the old object.  This
** means you can have several objects all pointing at the same thing.  When
** you start destroying objects, you don't want to destroy the original thing
** until you only have one reference left.
** ++*)
*/

  matrix::matrix(matrix& x)
{
  x.p->n++;    /* we're adding another reference, so increase */
  p = x.p;    /* the count. */
}


/*
-*++ ch_find(): search an input stream for a char, ignoring comments
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed: searches an input stream for c, ignoring comments starting
** with '#' to end of line.  Errfunc is the address of an error function,
** fname is  the name of the source file.
** ++*)
*/

    void ch_find(istream& from, char ch, void (*errfunc)(char *), char * fname)
{
  char c;
  UPCASE(ch);
  while(TRUE){
    while(from >> c, UPCASE(c), c != ch && c != '#')
      if(from.eof()) errfunc(fname); /* find ch or comment start */
    if (c == ch) break;
    /* If we made it here, a '#' comment start was found */
    while (from.get(c), c != EOL) /* 'get' doesn't skip whitespace */
      if(from.eof()) errfunc(fname); /* find end-of-line */
  }
}


/*
-*++ matrix::matrix(): read from a "standard" matrix file
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

  matrix::matrix(char * initfile) /* read from "standard" matrix file */
{
  void ch_find(istream& from, char ch, void (*errfunc)(char *), char * fname);
  void nonstandard(char *);
  filebuf f1;
  p = new matrep;
  if (f1.open(initfile,input) == 0)
    error2("cannot open matrix initializer file",initfile);
  istream from(&f1);

    /* Parse file initialization header  */

  ch_find(from,'r',nonstandard,initfile); /* find the 'r'... */
  ch_find(from,'=',nonstandard, initfile); /* ...followed by an `=` */
  from >> p->m_rows;        /* get row size */
  ch_find(from,'c',nonstandard, initfile); /* find the 'p'... */
  ch_find(from,'=',nonstandard, initfile); /* ...followed by an '=' */
  from >> p->m_cols;        /* get column size */
  ch_find(from,':',nonstandard, initfile); /* data follows ::: */
  ch_find(from,':',nonstandard, initfile);
  ch_find(from,':',nonstandard, initfile);

  p->m_vec = new _Vector*[rows()];
  for (int i = 0; i < rows(); i++)
    p->m_vec[i] = new _Vector(cols());
  p->n = 1;            /* we just created it: only 1 reference */
  for (int row = 0; row < rows(); row++)
    for(int col = 0; col < cols(); col++){
      from >> (*this)[row][col];
      if(from.bad()) error2("problem with matrix initializer file",initfile);
    }
}


/*
-*++ nonstandard(): an error message when a "non-standard" matrix file is found
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

void nonstandard(char * fname){
  cout << fname << " is a `non-standard' file.\n";
  cout << "A `standard' matrix file must start with\n";
  cout << "the dimensions of the matrix, i.e.:\n";
  cout << "rows=12 columns=14\n";
  cout << "or abbreviated:\n";
  cout << "r=12 c=14\n";
  cout << "comments follow `#' signs to end of line\n";
  cout << "data follows :::\n";
  exit(1);
}


/*
-*++ matrix::matrix(): create a special square matrix
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed: create a special square matrix.  If the "flag" string starts
** with:
**   I: identity matrix
**   R: square matrix with random values
** ++*)
*/

  matrix::matrix(char * flag, int dimension)
{
  double drand48();
  int i,j;
  if (flag[0] != 'I' && flag[0] != 'R')
    error("to create identity: \"I\"; random: \"R\"");
  p = new matrep;
  p->m_vec = new _Vector*[dimension];
  p->m_rows = dimension;
  p->m_cols = dimension;
  for ( i = 0; i < dimension; i++)
    p->m_vec[i] = new _Vector(dimension);
  p->n = 1;            /* so far, only one reference to this data */

  switch(flag[0]) {
  case 'I':
    for (i=0; i < rows(); i++)
      for ( j = 0; j < cols(); j++)
    (*this)[i][j] = (i == j ? 1 : 0);
    break;
  case 'R':
    for ( i=0; i < rows(); i++)
      for ( j = 0; j < cols(); j++)
    (*this)[i][j] = drand48();
  }
}


/*
-*++ matrix::~matrix(): matrix destructor.  Probably needs examining
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

 matrix::~matrix()   // this is probably wrong.
{
  if (--p->n == 0) {
    delete p->m_vec;        /* delete data */
    delete p;
  }
}


/*
-*++ matrix::write_standard(): write to a file in "standard" format
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed: write this matrix to filename in the 'standard' format
** ++*)
*/

void matrix::write_standard(char * filename, char * msg)
{
#ifndef DOS
  char *getwd(char *pathname);
#else DOS
  char *getcwd(char *pathname);
#endif
  char pathname[100];        /* for working directory information */
  filebuf f1;
  if (f1.open(filename,output) == 0)
    error2("cannot open or create matrix output file",filename);
  ostream to(&f1);
  to << "# " << filename << ": matrix file " << " written in `standard' format\n";
  long clock = time(0);
  to << "# " << asctime(localtime(&clock));
#ifndef DOS
  to << "# working directory: " << getwd(pathname) << "\n";
#else DOS
  to << "# working directory: " << getcwd(pathname) << "\n";
#endif
  to << "# " << msg << "\n";
  to << "\nRows=" << rows() << "  Columns=" << cols() << "\n";
  to << ":::\n";
  for (int row = 0; row < rows(); row++){
    for(int col = 0; col < cols(); col++){
      to << form("%6.6g  ",(*this)[row][col]);
      if(to.bad()) error2("problem with matrix output file",filename);
    }
    to << "\n";
  }
}


/*
-*++ matrix::operator=(): matrix "equals" using references
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

matrix &  matrix::operator=(matrix& rval)
{
  rval.p->n++;            /* tell the rval it has another reference to it */
  if(--p->n == 0) {        /* If nobody else is referencing us...*/
    delete p->m_vec;        /* ...nobody else knows to make us go away...*/
    delete p;
  }
  p = rval.p;     /* ...and we're pointing somewhere else (reference into rval) */
  return *this;
}


/*
-*++ matrix::operator+(): add two matrices
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

matrix & matrix::operator+(matrix& arg)
{
  if((rows() != arg.rows()) || (cols() != arg.cols()))
    error("matrices must be of the same dimension for addition!");
  matrix & sum= *new matrix(rows(),cols());
  for (int i=0; i< rows(); i++)
    for (int j = 0; j < cols(); j++)
      sum[i][j] = (*this)[i][j] + arg[i][j];
  return sum;
}


/*
-*++ matrix::operator+(): add a double value to each element
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

matrix & matrix::operator+(double arg)
{
  matrix & sum= *new matrix(rows(),cols());
  for (int i=0; i< rows(); i++)
    for (int j = 0; j < cols(); j++)
      sum[i][j] = (*this)[i][j] + arg;
  return sum;
}


/*
-*++ matrix::operator*(): multiply each element by a double value
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

matrix & matrix::operator*(double arg)
{
  matrix & sum= *new matrix(rows(),cols());
  for (int i=0; i< rows(); i++)
    for (int j = 0; j < cols(); j++)
      sum[i][j] = (*this)[i][j] * arg;
  return sum;
}


/*
-*++ matrix::operator*(): matrix multiply
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
**      9 Dec 87    Bruce Eckel    Changed the test -- I had arg1.rows !=
** arg2.cols, which was backwards
** ++*)
**
** (*++ detailed:
** ++*)
*/

matrix & matrix::operator*(matrix& arg)
{
  if(cols() != arg.rows())
    error("# rows of second mat must equal # cols of first for multiply!");
  matrix & result= *new matrix(rows(),arg.cols());
  for(int row = 0; row < rows(); row++)
    for(int col = 0; col < arg.cols(); col++){
      double sum = 0;
      for(int i = 0; i < cols(); i++)
      sum += (*this)[row][i] * arg[i][col];
      result[row][col] = sum;
    };
  return result;
}


/*
-*++ matrix::operator-(): subtract one matrix from another
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

matrix & matrix::operator-(matrix& arg)
{
    if((rows() != arg.rows()) || (cols() != arg.cols()))
    error("matrices must be of the same dimension for subtraction!\n");
    matrix & sum= *new matrix(rows(),cols());
    for (int i=0; i< rows(); i++)
    for (int j = 0; j < cols(); j++)
        sum[i][j] = (*this)[i][j] - arg[i][j];
    return sum;
}


/*
-*++ matrix::operator-(): unary minus each element in a matrix
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

matrix & matrix::operator-()
{
    for (int i=0; i< rows(); i++)
    for (int j = 0; j < cols(); j++)
        (*this)[i][j] *= -1;
    return *this;
}


/*
-*++ matrix::operator==(): test for matrix eqivalence
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

int matrix::operator==(matrix& arg)
{
  if((rows() != arg.rows()) || (cols() != arg.cols()))
    error("matrices must be of the same dimension for equivalence test!\n");
  for (int i=0; i< rows(); i++)
    for (int j = 0; j < cols(); j++)
      if( (*this)[i][j] != arg[i][j])
    return 0;
  return ~0;
}


/*
-*++ matrix::operator!=(): test for matrix eqivalence
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

int matrix::operator!=(matrix& arg)
{
  if((rows() != arg.rows()) || (cols() != arg.cols()))
    error("matrices must be of the same dimension for inequivalence!\n");
  for (int i=0; i< rows(); i++)
    for (int j = 0; j < cols(); j++)
      if( (*this)[i][j] == arg[i][j])
    return 0;
  return ~0;
}


/*
-*++ matrix::operator<<: print a matrix to standard out
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

ostream& operator<<(ostream &s, matrix& mat)
{
  for (int i=0; i< mat.rows(); i++){
    for (int j = 0; j < mat.cols(); j++)
      s <<  form("%6.6g  ",mat[i][j]);
    s << "\n";
  }
  return s;
}


/*
-*++ matrix::operator>>():
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

istream& operator>>(istream& s, matrix& mat)
{
  double val = 0;
  for (int row = 0; row < mat.rows(); row++)
    for(int col = 0; col < mat.cols(); col++){
      s >> val;
      if (s.bad())
    mat.error("bad matrix initialization input\n");
      if (s.eof())
    mat.error("end of matrix initialization input before matrix full\n");
      mat[row][col] = val;
    }
  return s;
}


/*
-*++ matrix::file_get(): fills a matrix from a `non-standard' ASCII file
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed: fills a matrix of known size from a `non-standard' file of
** numbers
** ++*)
*/

void matrix::file_get(char * infile)
{
  filebuf f1;
  if (f1.open(infile,input) == 0)
    error2("cannot open matrix initializer file",infile);
  istream from(&f1);
  for (int row = 0; row < rows(); row++)
    for(int col = 0; col < cols(); col++){
      from >> (*this)[row][col];
      if(from.bad()) error2("problem with matrix initializer file",infile);
    }
}


/*
-*++ matrix::file_put(): writes a matrix to a "non-standard" ASCII file
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

void matrix::file_put(char * outfile)
{
  filebuf f1;
  if (f1.open(outfile,output) == 0)
    error2("cannot open or create matrix output file",outfile);
  ostream to(&f1);
  for (int row = 0; row < rows(); row++){
    for(int col = 0; col < cols(); col++){
      to << (*this)[row][col];
      to << "  ";
      if(to.bad()) error2("problem with matrix output file",outfile);
    }
    to << "\n";
  }
}


/*
-*++ operator<<(): output for XY_vec
**
** (*++ history:
**     11 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

ostream& operator<<(ostream &s, XY_vec& v) {
    s << " X:       Y:\n";
    for(int i = 0; i < v.size(); i++)
    s << form("%6.6g   %6.6g\n",v.X.Indep(i),v.Y[i]);
    return s;
}


/*
-*++ XY_vec::floatvect(): return a pointer to an array of floating-point #s
**
** (*++ history:
**     14 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed: What is this terminated by?
** ++*)
*/

float *  XY_vec::floatvect()
{
    float * result = new float[2 * size()];
    for (int i = 0; i < size() * 2; i += 2) {
    result[i] = (float)X.Indep(i/2);
    result[i+1] = (float)Y[i/2];
    }
    return result;
}


/*
-*++ ColVec::doublvect(): returns a simple array of double numbers
**
** (*++ history:
**     15 Jan 88    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

double * ColVec::doublvect()
{
    double * result = new double[size()];
    for (int i = 0; i < size(); i++)
    result[i] = (*this)[i];
    return result;
}